#include "AnaDef.h"
#include "XGLUtils.h"

int main(int argc, char * argv[])
{
  for(Int_t ii=0; ii<argc; ii++){
    printf("%d: %s\n",ii, argv[ii]);
  }

  if(argc!=2){
    printf("argc!=2\n");return 1;
  }

  const Int_t kmode=atoi(argv[1]);

  TList *ll=new TList;
  TGraph * grtrueN=new TGraph;
  grtrueN->SetName("trueN");
  ll->Add(grtrueN);

  TGraph * grfitN=new TGraph;
  grfitN->SetName("fitN");
  ll->Add(grfitN);

  TGraphErrors * grratioN=new TGraphErrors;
  grratioN->SetName("ratioN");
  ll->Add(grratioN);
  
  TRandom3 ran(1);
  Int_t ifit = 0;
  for(Int_t nn = 100000; nn> 100; nn/=2){
    printf("nn %d\n",nn);

    TH1D * hh=new TH1D(Form("hh%d",ifit),"",100,100,200);
    ll->Add(hh);
    for(Int_t ii=0; ii<nn; ii++){
      hh->Fill(ran.Gaus(150,5));
    }

    grtrueN->SetPoint(ifit, nn, nn);

    if(kmode==0){
      TF1 f1("f1","gausn",100,200);
      f1.SetParameters(1,hh->GetMean(), hh->GetRMS());
      hh->Fit(&f1);//,"L");
      grfitN->SetPoint(ifit,nn,f1.GetParameter(0));
      grratioN->SetPoint(ifit,nn,f1.GetParameter(0)/nn);
      const Double_t err = f1.GetParError(0);
      grratioN->SetPointError(ifit,0,err/nn);
    }
    else if(kmode==1){
      TF1 f1("f1","gausn",100,200);
      f1.SetParameters(1,hh->GetMean(), hh->GetRMS());
      hh->Fit(&f1,"L");
      grfitN->SetPoint(ifit,nn,f1.GetParameter(0));
      grratioN->SetPoint(ifit,nn,f1.GetParameter(0)/nn);
      const Double_t err = f1.GetParError(0);
      grratioN->SetPointError(ifit,0,err/nn);
    }
    else if(kmode==2){
      Double_t pars[]={1,hh->GetMean(), hh->GetRMS()};
      Double_t errs[10];
      XGLUtils::FFunc func = XGLUtils::RawGaus;
      const Int_t npar = 3;
      const Int_t kfail = XGLUtils::ChisquareFit(hh, func, npar, pars, errs);
      printf("kmode %d done with kfail %d\n\n", kmode, kfail);
      grfitN->SetPoint(ifit,nn,pars[0]);
      grratioN->SetPoint(ifit,nn,pars[0]/nn);
      grratioN->SetPointError(ifit,0, errs[0]/nn);

      TH1D *hfit=XGLUtils::GetHfit(Form("%sfit",hh->GetName()), func, pars, hh->GetXaxis()->GetXmin(), hh->GetXaxis()->GetXmax());
      hfit->SetLineColor(kRed);
      ll->Add(hfit);
    }
    else if(kmode==3){
      //un-normalized hist -- binned likelihood
      Double_t pars[]={1, hh->GetMean(), hh->GetRMS()};
      Double_t errs[10];
      XGLUtils::FFunc func = XGLUtils::RawGaus;
      const Int_t npar = 3;
      const Int_t kfail = XGLUtils::BinnedLikelihoodFit(hh, func, npar, pars, errs);
      printf("kmode %d done with kfail %d\n\n", kmode, kfail);
      grfitN->SetPoint(ifit,nn,pars[0]);
      grratioN->SetPoint(ifit,nn,pars[0]/nn);
      grratioN->SetPointError(ifit,0, errs[0]/nn);

      TH1D *hfit=XGLUtils::GetHfit(Form("%sfit",hh->GetName()), func, pars, hh->GetXaxis()->GetXmin(), hh->GetXaxis()->GetXmax());
      hfit->SetLineColor(kRed);
      ll->Add(hfit);
    }
    else if(kmode==4){
      //normalized hist -- chi2
      hh->Sumw2();//important!!
      hh->Scale(1./hh->Integral(0,100000));
      Double_t pars[]={1,hh->GetMean(), hh->GetRMS()};
      Double_t errs[10];
      XGLUtils::FFunc func = XGLUtils::RawGaus;
      const Int_t npar = 3;
      const Int_t kfail = XGLUtils::ChisquareFit(hh, func, npar, pars, errs);
      printf("kmode %d done with kfail %d\n\n", kmode, kfail);
      grfitN->SetPoint(ifit,nn,pars[0]*nn);
      grratioN->SetPoint(ifit,nn,pars[0]);
      grratioN->SetPointError(ifit,0, errs[0]);

      TH1D *hfit=XGLUtils::GetHfit(Form("%sfit",hh->GetName()), func, pars, hh->GetXaxis()->GetXmin(), hh->GetXaxis()->GetXmax());
      hfit->SetLineColor(kRed);
      ll->Add(hfit);
    }
    else if(kmode==5){
      //normalized hist -- binned likelihood
      hh->Sumw2();//important!!
      hh->Scale(1./hh->Integral(0,100000));
      Double_t pars[]={1./hh->Integral(0,100000), hh->GetMean(), hh->GetRMS()};
      Double_t errs[10];
      XGLUtils::FFunc func = XGLUtils::RawGaus;
      const Int_t npar = 3;
      const Int_t kfail = XGLUtils::BinnedLikelihoodFit(hh, func, npar, pars, errs);
      printf("kmode %d done with kfail %d\n\n", kmode, kfail);

      grfitN->SetPoint(ifit,nn,pars[0]*nn);
      grratioN->SetPoint(ifit,nn,pars[0]);
      grratioN->SetPointError(ifit,0, errs[0]);

      TH1D *hfit=XGLUtils::GetHfit(Form("%sfit",hh->GetName()), func, pars, hh->GetXaxis()->GetXmin(), hh->GetXaxis()->GetXmax());
      hfit->SetLineColor(kRed);
      ll->Add(hfit);
    }
    else{
      printf("kmode wrong!! %d\n", kmode);exit(1);
    }

    printf("fit result %d %d %e %e %e\n", ifit, nn, grfitN->GetY()[ifit], grratioN->GetY()[ifit], grratioN->GetEY()[ifit]);
    ifit++;
  }

  grratioN->SetMarkerStyle(20);
  grratioN->Draw("ap");
  gPad->SetLogx();

  TFile *fout=new TFile(Form("outDemoFit_mode%d.root", kmode),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();

  return 0;
}
